load packages

library(tidyverse)
library(knitr)
devtools::install_github("dcosme/specr", ref = "plotmods")
library(specr)

define palettes

pal_outcome = wesanderson::wes_palette("Zissou1", 3, "continuous")
pal_sca = c(pal_outcome[3], pal_outcome[3], "grey")
pal_condition = wesanderson::wes_palette("Zissou1", 16, "continuous")
pal_condition = c(pal_condition[1], pal_condition[3],
                  pal_condition[14], pal_condition[5],
                  pal_condition[15], pal_condition[16],
                  pal_condition[9], pal_condition[12])

set seed

set.seed(6523)

define SCA wrapper functions

plot functions

plot_sca = function(data, combined = TRUE, labels = c("A", "B"), title = FALSE, limits = NULL,
                    point_size = 1, point_alpha = 1, ci_alpha = .5, ci_size = .5, 
                    text_size = 14, title_size = 6,
                    choices  = c("x", "y", "test", "stimulus", "contrast", "controls"),
                    remove_y = FALSE, remove_facet = FALSE) {
  
  if (combined == TRUE) {
      p1 = plot_curve(data, point_size = point_size, point_alpha = point_alpha,
                      ci_alpha = ci_alpha, ci_size = ci_size,
                      limits = limits) +
        geom_hline(aes(yintercept = median(filter(data, y == "body fat percentage")$estimate), linetype = "body fat percentage"), color = pal_outcome[1], size = 1) +
        geom_hline(aes(yintercept = median(filter(data, y == "BMI")$estimate), linetype = "BMI"), color = pal_outcome[2], size = 1) +
        geom_hline(aes(yintercept = median(filter(data, y == "enactment")$estimate), linetype = "enactment"), color = pal_outcome[3], size = 1) +
        scale_linetype_manual(name = "", values = c(2, 2, 2), guide = guide_legend(override.aes = list(color = c(pal_outcome[1], pal_outcome[2], pal_outcome[3])))) +
        labs(x = "", y = "standarized\nregression coefficient\n")  +
        theme(legend.position = "top",
              text = element_text(size = text_size))
      
      if (title == TRUE) {
        if (is.null(limits)) {
          p1 = p1 + annotate("text", -Inf, Inf, label = unique(data$x), fontface = 2, size = title_size,
                       x = 0.5*(1 + nrow(data)), 
                       y = max(data$conf.high))
        } else {
          p1 = p1 + annotate("text", -Inf, Inf, label = unique(data$x), fontface = 2, size = title_size,
                       x = 0.5*(1 + nrow(data)), 
                       y = limits[2])
        }
      }
      
      p2 = plot_choices(data, choices = choices,
                        ignore_vars = c("ROI", "neural signature"), rename_controls = "covariates") +
        labs(x = "\nspecifications (ranked)")  +
        theme(strip.text.x = element_blank(),
              text = element_text(size = text_size))
        
  }
  else {
      p1 = plot_curve(data, point_size = point_size, point_alpha = point_alpha,
                      ci_alpha = ci_alpha, ci_size = ci_size) +
        geom_hline(yintercept = 0, linetype = "solid", color = "black", size = .5) +
        labs(x = "", y = "standarized\nregression coefficient\n") +
        theme(text = element_text(size = text_size))
      
      if (title == TRUE) {
        if (is.null(limits)) {
          p1 = p1 + annotate("text", -Inf, Inf, label = unique(data$y), fontface = 2, size = title_size,
                       x = 0.5*(1 + nrow(data)), 
                       y = max(data$conf.high))
        } else {
          p1 = p1 + annotate("text", -Inf, Inf, label = unique(data$y), fontface = 2, size = title_size,
                       x = 0.5*(1 + nrow(data)), 
                       y = limits[2])
        }
      }
      
      p2 = plot_choices(data, choices = choices,
                        ignore_vars = c("ROI", "neural signature", "not included"), rename_controls = "covariates") +
        labs(x = "\nspecification number (ranked)") +
        theme(strip.text.x = element_blank(),
              text = element_text(size = text_size))
  }
  
  if (remove_y == TRUE) {
    p1 = p1 + labs(y = "")
    
    p2 = p2 + theme(axis.text.y = element_blank(),
                    axis.ticks.y = element_blank()) +
      labs(y = "")
  }

  if (remove_facet == TRUE) {
    p2 = p2 + theme(strip.text.y = element_blank())
    }
  
  plot_specs(plot_a = p1,
             plot_b = p2,
             labels = labels,
             rel_height = c(1, 2))
}

plot_sca_compare = function(data, pointrange = TRUE, labels = c("A", "B"), 
                            rel_heights = c(.75, .25), rel_widths = c(.75, .25), 
                            title = FALSE, text_size = 14, title_size = 6, n_rows = 1,
                            remove_x = FALSE, remove_y = FALSE, sig = NULL) {
  
  # define median bootstrapping function
  median_cl_boot = function(x, conf = 0.95, df = TRUE, ci = "low") {
  
    lconf = (1 - conf)/2
    uconf = 1 - lconf
    require(boot)
    bmedian = function(x, ind) median(x[ind])
    bt = boot(x, bmedian, 1000)
    bb = boot.ci(bt, type = "perc")
    
    if (df == TRUE){
      data.frame(y = median(x),
                 ymin = quantile(bt$t, lconf), 
                 ymax = quantile(bt$t, uconf))
      
    } else {
      if (ci == "low") {
        quantile(bt$t, lconf)
      } else {
        quantile(bt$t, uconf)
      }
    }
  }
    
  # merge and tidy for plotting
  plot.data = data %>%
    group_by(x) %>%
    arrange(estimate) %>%
    mutate(specification = row_number()) %>%
    ungroup() %>%
    unique()
  
  # labels
  labs = plot.data %>%
    group_by(x) %>%
    summarize(med = median(estimate),
              low = median_cl_boot(estimate, df = FALSE, ci = "low"),
              high = median_cl_boot(estimate, df = FALSE, ci = "high")) %>%
    mutate(range = max(high) - min(low),
           estimate = ifelse(med > 0, high + (range / 10), low - (range / 10)),
           label = ifelse(x %in% sig, "*", ""))
  
  # plot curves
  if (pointrange == TRUE) {
    a = plot.data %>%
    ggplot(aes(specification, estimate, color = x)) +
      geom_linerange(aes(ymin = conf.low, ymax = conf.high), size = .1) +
      geom_point() +
      geom_hline(yintercept = 0, linetype = "solid", color = "black", size = 1) +
      scale_color_manual(name = "", values = pal_condition) +
      labs(x = "\nspecification number (ranked)", y = "standarized\negression coefficient\n") + 
      theme_minimal() + 
      theme(strip.text = element_blank(), 
            axis.line = element_line("black", size = 0.5), 
            legend.position = c(.5, .1), 
            legend.direction = "horizontal",
            panel.spacing = unit(0.75, "lines"), 
            axis.text = element_text(colour = "black"),
            text = element_text(size = text_size))
    if (title == TRUE) {
      a = a + annotate("text", -Inf, Inf, label = unique(plot.data$y), fontface = 2, size = title_size,
                       x = 0.5*(min(plot.data$specification) + max(plot.data$specification)), 
                       y = max(plot.data$conf.high))
    }
    
  } else {
    a = plot.data %>%
      ggplot(aes(specification, estimate, color = x)) +
      geom_point() +
      geom_hline(yintercept = 0, linetype = "solid", color = "black", size = 1) +
      scale_color_manual(name = "", values = pal_condition) +
      labs(x = "\nspecification number (ranked)", y = "standarized\nregression coefficient\n") + 
      theme_minimal() + 
      theme(strip.text = element_blank(), 
            axis.line = element_line("black", size = 0.5), 
            legend.position = "none", 
            legend.direction = "horizontal",
            panel.spacing = unit(0.75, "lines"), 
            axis.text = element_text(colour = "black"),
            text = element_text(size = text_size))
    if (title == TRUE) {
      a = a + annotate("text", -Inf, Inf, label = unique(plot.data$y), fontface = 2, size = title_size,
                       x = 0.5*(min(plot.data$specification) + max(plot.data$specification)), 
                       y = max(plot.data$estimate))    
      }
  }
  
    b = plot.data %>%
      group_by(x) %>%
      mutate(order = median(estimate)) %>%
      ggplot(aes(reorder(x, order), estimate, fill = x)) +
      stat_summary(fun.y = "median", geom = "bar") +
      stat_summary(fun.data = median_cl_boot, geom = "errorbar", width = 0) +
      geom_text(data = labs, aes(label = label, x = x, y = estimate), size = 6) +
      scale_fill_manual(name = "", values = pal_condition) +
      scale_y_continuous(breaks = scales::pretty_breaks(n = 4)) + 
      labs(x = "\n", y = "median effect\n") + 
      theme_minimal() + 
      theme(strip.text = element_blank(), 
            axis.line = element_line("black", size = 0.5), 
            legend.position = "none", 
            panel.spacing = unit(0.75, "lines"), 
            axis.text = element_text(colour = "black"),
            text = element_text(size = text_size),
            axis.text.x = element_text(angle = 45, hjust = 1))
    
  if (n_rows == 1) {
    a = a + theme(legend.position = c(.5, .1))
    b = b + coord_flip() +
      labs(x = "\n", y = "\nmedian effect") + 
      theme(axis.text.x = element_text(angle = 0, hjust = 1),
            axis.text.y = element_blank())
  }     
    

  if (remove_x == TRUE) {
    a = a + labs(x = "")
    
    if (n_rows == 1) {
      b = b + labs(y = "")
    } else {
      b = b + labs(x = "")
    }
  }    
  
  if (remove_y == TRUE) {
    a = a + labs(y = "")
    
    if (n_rows == 1) {
      b = b + labs(x = "")
    } else {
      b = b + labs(y = "")
    }
  }  
    
  cowplot::plot_grid(a, b, labels = labels, rel_heights = rel_heights, rel_widths = rel_widths, nrow = n_rows)
}

summary and bootstapping functions

Process overview

  • shuffle value-std before spreading
  • run SCA with shuffled values
  • extract median estimate, % positive (& significant p < .05), % negative (& signiticant p < .05)
  • repeat process 500 times
run_boot = function(data, var, n_shuffles, rerun) {
  
  if (rerun == FALSE & file.exists(sprintf("boot/boot_%s.RDS", var))) {
    out = readRDS(sprintf("boot/boot_%s.RDS", var))
  } else {
    out = sca_boot(data = data_combined, var = var, n_shuffles = n_shuffles)
    saveRDS(out, sprintf("boot/boot_%s.RDS", var))
  }
  return(out)
  
}

sca_boot = function(data, var, n_shuffles) {
  
  # set seed
  set.seed(63)
  
  # initialize data frame
  df = data.frame()
  
  for (i in 1:n_shuffles){
    # shuffle `value-std`
    data_shuffled = data %>%
      mutate(shuffled = sample(value_std)) %>%
      select(-value_std)
    
    # separate uniformity and association data
    data1_shuffled = data_shuffled %>%
      filter(grepl("uniformity|ROI", test)) %>%
      select(-test) %>%
      spread(process, shuffled)
    
    data2_shuffled = data_shuffled %>%
      filter(grepl("association|ROI", test)) %>%
      select(-test) %>%
      spread(process, shuffled)
    
    # define control variables
    controls = unique(data_shuffled$process)[unique(data_shuffled$process) != var]
      
    # run scas
    results_data1 = run_specs(df = data1_shuffled,
                              y = c("bmi", "fat", "enact_prop"), 
                              x = var,
                              controls = controls,
                              model = c("lm"), 
                              subsets = list(stimulus = unique(data1_shuffled$stimulus),
                                             contrast = unique(data1_shuffled$contrast))) %>%
      mutate(test = ifelse(grepl("ROI", controls), "ROI", 
                    ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI", 
                    ifelse(controls == "craving_regulation_PEV", "neural signature",     
                    ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "uniformity")))))
    
    results_data2 = run_specs(df = data2_shuffled,
                              y = c("bmi", "fat", "enact_prop"), 
                              x = var,
                              controls = controls,
                              model = c("lm"), 
                              subsets = list(stimulus = unique(data2_shuffled$stimulus),
                                             contrast = unique(data2_shuffled$contrast))) %>%
      mutate(test = ifelse(grepl("ROI", controls), "ROI", 
                    ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI", 
                    ifelse(controls == "craving_regulation_PEV", "neural signature",     
                    ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "association")))))
    
    # merge and extract info
    results = bind_rows(results_data1, results_data2) %>%
      mutate(duplicated = ifelse(duplicated(estimate) & duplicated(std.error), 1, 0)) %>%
      filter(!duplicated) %>%
      unique() %>%
      ungroup() %>%
      mutate(pos = ifelse(estimate > 0, 1, 0),
             neg = ifelse(estimate < 0, 1, 0),
             sig = ifelse(p.value < .05, 1, 0),
             pos_sig = ifelse(pos == 1 & sig == 1, 1, 0),
             neg_sig = ifelse(neg == 1 & sig == 1, 1, 0)) %>%
      group_by(x, y) %>%
      summarize(median = median(estimate),
                n = n(),
                n_positive = sum(pos),
                n_negative = sum(neg),
                n_significant = sum(sig),
                n_positive_sig = sum(pos_sig),
                n_negative_sig = sum(neg_sig)) %>%
      mutate(shuffle = i) %>%
      select(shuffle, everything()) %>%
      ungroup() %>%
      mutate_if(is.character, stringr::str_replace_all, pattern = "enact_prop", replacement = "enactment") %>%
      mutate_if(is.character, stringr::str_replace_all, pattern = "bmi", replacement = "BMI") %>%
      mutate_if(is.character, stringr::str_replace_all, pattern = "fat", replacement = "body fat percentage") %>%
      mutate_if(is.character, stringr::str_replace_all, pattern = "_", replacement = " ")
  
    # join to data frame
    df = rbind(df, as.data.frame(results))
    rm(results)
    
  }
  
  return(df)
  
}

run_sca = function(data, var, outcome, boot = FALSE, n_shuffles = 1, rerun = FALSE) {
  
    # define median bootstrapping function
    median_cl_boot = function(estimate, conf = 0.95) {
  
      lconf = (1 - conf)/2
      uconf = 1 - lconf
      require(boot)
      bmedian = function(estimate, ind) median(estimate[ind])
      bt = boot(estimate, bmedian, 1000)
      bb = boot.ci(bt, type = "perc")
      data.frame(obs_median = median(estimate), 
                 obs_conf_low = quantile(bt$t, lconf), 
                 obs_conf_high = quantile(bt$t, uconf))
    
    }
    
    # separate uniformity and association data
    data1 = data %>%
      filter(grepl("uniformity|ROI", test)) %>%
      select(-test) %>%
      spread(process, value_std)
    
    data2 = data %>%
      filter(grepl("association|ROI", test)) %>%
      select(-test) %>%
      spread(process, value_std)
    
    # define control variables
    if (length(var) > 1) {
      controls = unique(data$process)
    } else {
      controls = unique(data$process)[!unique(data$process) %in% var]
    }
      
    # run scas
    results_data1 = run_specs(df = data1,
                              y = outcome, 
                              x = var,
                              controls = controls,
                              model = c("lm"), 
                              subsets = list(stimulus = unique(data1$stimulus),
                                             contrast = unique(data1$contrast))) %>%
      mutate(test = ifelse(grepl("ROI", controls), "ROI", 
                    ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI", 
                    ifelse(controls == "craving_regulation_PEV", "neural signature",     
                    ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "uniformity")))))
    
    results_data2 = run_specs(df = data2,
                              y = outcome, 
                              x = var,
                              controls = controls,
                              model = c("lm"), 
                              subsets = list(stimulus = unique(data2$stimulus),
                                             contrast = unique(data2$contrast)))  %>%
      mutate(test = ifelse(grepl("ROI", controls), "ROI", 
                    ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI", 
                    ifelse(controls == "craving_regulation_PEV", "neural signature",     
                    ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "association")))))
    
    # merge and extract info & exclude duplicates
    results = bind_rows(results_data1, results_data2) %>%
      filter(!x == controls) %>%
      mutate(stimulus = ifelse(grepl("snack", subsets), "snack",
                        ifelse(grepl("meal", subsets), "meal",
                        ifelse(grepl("dessert", subsets), "dessert",
                        ifelse(grepl("food", subsets), "food", "all")))),
             contrast = ifelse(grepl("nature", subsets), "nature",
                        ifelse(grepl("rest", subsets), "rest", "all")),
             duplicated = ifelse(duplicated(estimate) & duplicated(std.error), 1, 0)) %>%
      filter(!duplicated) %>%
      unique() %>%
      ungroup() %>%
      mutate_if(is.character, stringr::str_replace_all, pattern = "enact_prop", replacement = "enactment") %>%
      mutate_if(is.character, stringr::str_replace_all, pattern = "bmi", replacement = "BMI") %>%
      mutate_if(is.character, stringr::str_replace_all, pattern = "fat", replacement = "body fat percentage") %>%
      mutate_if(is.character, stringr::str_replace_all, pattern = "_", replacement = " ")
    
    median_ci = results %>%
      group_by(x, y) %>%
      do({
        median_cl_boot(.$estimate)
      })
        
    summary = results %>%
      mutate(pos = ifelse(estimate > 0, 1, 0),
             neg = ifelse(estimate < 0, 1, 0),
             sig = ifelse(p.value < .05, 1, 0),
             pos_sig = ifelse(pos == 1 & sig == 1, 1, 0),
             neg_sig = ifelse(neg == 1 & sig == 1, 1, 0)) %>%
      group_by(x, y) %>%
      summarize(obs_n = n(),
                obs_n_positive = sum(pos),
                obs_n_negative = sum(neg),
                obs_n_significant = sum(sig),
                obs_n_positive_sig = sum(pos_sig),
                obs_n_negative_sig = sum(neg_sig)) %>%
      left_join(., median_ci, by = c("x", "y")) %>%
      select(x, y, obs_median, obs_conf_low, obs_conf_high, everything())
      

  
    if (boot == TRUE) {
      boot_out = run_boot(data, var, n_shuffles, rerun)
      return(list(results = results, summary = summary, boot = boot_out))
    } else {
      return(list(results = results, summary = summary))
    }
    
}

boot_stats = function(data) {
  
  data$boot %>%
    mutate(y = ifelse(grepl("%", y), "body fat percentage", y)) %>%
    left_join(., data$summary, c("x", "y")) %>%
    mutate(extreme_median = ifelse(obs_median > 0 & median >= obs_median, 1,
                            ifelse(obs_median < 0 & median <= obs_median, 1, 0)),
           extreme_positive = ifelse(n_positive >= obs_n_positive, 1, 0),
           extreme_negative = ifelse(n_negative >= obs_n_negative, 1, 0),
           extreme_positive_sig = ifelse(n_positive_sig >= obs_n_positive_sig, 1, 0),
           extreme_negative_sig = ifelse(n_negative_sig >= obs_n_negative_sig, 1, 0)) %>%
    group_by(x, y) %>%
    summarize(n = n(),
              extreme_median = sum(extreme_median),
              extreme_positive = sum(extreme_positive),
              extreme_negative = sum(extreme_negative),
              extreme_positive_sig = sum(extreme_positive_sig),
              extreme_negative_sig = sum(extreme_negative_sig)) %>%
    gather(variable, n_extreme, contains("extreme")) %>%
    mutate(p_value = n_extreme / n,
           p_value = ifelse(p_value == 1.000, "1.000", 
               ifelse(p_value < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p_value))))) %>%
    select(x, y, variable, p_value) %>%
    spread(variable, p_value) %>%
    left_join(., select(ungroup(data$summary), x, y, obs_median, obs_conf_low, obs_conf_high, obs_n_positive_sig, obs_n_negative_sig), by = c("x", "y")) %>%
    mutate(`Mdn [95% CI]` = sprintf("%.2f [%.2f, %.2f]", obs_median, obs_conf_low, obs_conf_high)) %>%
    select(y, `Mdn [95% CI]`, extreme_median, obs_n_positive_sig, extreme_positive_sig, obs_n_negative_sig, extreme_negative_sig) %>%
    rename("Mdn p" = extreme_median,
           "Positive N" = obs_n_positive_sig,
           "Positive p" = extreme_positive_sig,
           "Negative N" = obs_n_negative_sig,
           "Negative p" = extreme_negative_sig) %>%
    gather(var, val, -c(x, y)) %>%
    unite(y, y, var, sep = " ") %>%
    spread(y, val)
  
    # data$boot %>%
    # left_join(., data$summary, c("x", "y")) %>%
    # mutate(extreme_median = ifelse(obs_median > 0 & median >= obs_median, 1,
    #                         ifelse(obs_median < 0 & median <= obs_median, 1, 0)),
    #        extreme_positive = ifelse(n_positive >= obs_n_positive, 1, 0),
    #        extreme_negative = ifelse(n_negative >= obs_n_negative, 1, 0),
    #        extreme_positive_sig = ifelse(n_positive_sig >= obs_n_positive_sig, 1, 0),
    #        extreme_negative_sig = ifelse(n_negative_sig >= obs_n_negative_sig, 1, 0)) %>%
    # group_by(x, y) %>%
    # summarize(n = n(),
    #           extreme_median = sum(extreme_median),
    #           extreme_positive = sum(extreme_positive),
    #           extreme_negative = sum(extreme_negative),
    #           extreme_positive_sig = sum(extreme_positive_sig),
    #           extreme_negative_sig = sum(extreme_negative_sig)) %>%
    # gather(variable, n_extreme, contains("extreme")) %>%
    # mutate(p_value = n_extreme / n,
    #        se = sqrt((p_value * (1 - p_value)) / n),
    #        p_value = ifelse(p_value == 1.000, "1.000", 
    #            ifelse(p_value < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p_value))))) %>%
    # select(x, y, variable, p_value, se) %>%
    # left_join(., select(ungroup(data$summary), x, y, obs_median, obs_conf_low, obs_conf_high, obs_n_positive_sig, obs_n_negative_sig), by = c("x", "y")) %>%
    # mutate(`Mdn [95% CI]` = sprintf("%.2f [%.2f, %.2f]", obs_median, obs_conf_low, obs_conf_high)) %>%
    # #select(y, `Mdn [95% CI]`, extreme_median, obs_n_positive_sig, extreme_positive_sig, obs_n_negative_sig, extreme_negative_sig) %>%
    # mutate(variable = gsub("extreme_median", "Mdn p", variable),
    #        variable = gsub("obs_n_positive_sig", "Positive N", variable),
    #        variable = gsub("extreme_positive_sig", "Positive p", variable),
    #        variable = gsub("obs_n_negative_sig", "Negative N", variable),
    #        variable = gsub("extreme_negative_sig", "Negative p", variable),
    #        `Mdn [95% CI]2` = sprintf("%.2f [%.2f, %.2f]", obs_median, obs_median - (1.96*se), obs_median + (1.96*se))) %>%
    #   filter(!grepl("extreme", variable))
    # gather(var, val, -c(x, y)) %>%
    # unite(y, y, var, sep = " ") %>%
    # spread(y, val)
    # 
}

load and tidy data

Read betas_dataset.RDS and dots_dataset.RDS saved from the prediction_models.Rmd

  • Ran uniformity and association test data separately because the models were breaking when trying to include test type as a covariate
source("load_data.R")

ind_diffs1 = ind_diffs %>%
  mutate(bmi = as.numeric(scale(bmi)),
         fat = as.numeric(scale(fat)))

ema_enact1 = ema_enact %>%
  mutate(enact_prop = as.numeric(scale(enact_prop)))

betas_std = readRDS("betas_dataset.RDS") %>%
  filter(session == "all") %>%
  group_by(subjectID, process, condition, control) %>%
  mutate(value_std = mean(meanPE_std, na.rm = TRUE)) %>%
  ungroup() %>%
  select(subjectID, process, condition, control, value_std) %>%
  unique() %>%
  mutate(process = sprintf("%s_ROI", process)) %>%
  filter(grepl("snack|meal|dessert|food", condition))

data_unif = readRDS("dots_dataset.RDS") %>%
  filter(session == "all" & mask == "unmasked") %>%
  filter((test == "uniformity" & !process == "craving_regulation") | 
           (test == "association" & process == "craving_regulation"))  %>%
  mutate(process = sprintf("%s_PEV", process)) %>%
  rename("value_std" = dotProduct_std) %>%
  select(subjectID, process, condition, control, value_std) %>%
  unique() %>%
  filter(grepl("snack|meal|dessert|food", condition)) %>%
  bind_rows(betas_std) %>%
  left_join(., ind_diffs1, by = "subjectID") %>%
  left_join(., ema_enact1, by = "subjectID") %>%
  select(-c(sample, DBIC_ID, age, gender)) %>%
  rename("contrast" = control,
         "stimulus" = condition) %>%
  mutate(stimulus = as.factor(stimulus),
         contrast = as.factor(contrast)) %>%
  spread(process, value_std)

data_assoc = readRDS("dots_dataset.RDS") %>%
  filter(session == "all" & mask == "unmasked" & test == "association") %>%
  mutate(process = sprintf("%s_PEV", process)) %>%
  rename("value_std" = dotProduct_std) %>%
  select(subjectID, process, condition, control, value_std) %>%
  unique() %>%
  filter(grepl("snack|meal|dessert|food", condition)) %>%
  bind_rows(betas_std) %>%
  left_join(., ind_diffs1, by = "subjectID") %>%
  left_join(., ema_enact1, by = "subjectID") %>%
  select(-c(sample, DBIC_ID, age, gender)) %>%
  rename("contrast" = control,
         "stimulus" = condition) %>%
  mutate(stimulus = as.factor(stimulus),
         contrast = as.factor(contrast)) %>%
  spread(process, value_std)

data_unif_gathered = data_unif %>%
  gather(process, value_std, contains("ROI"), contains("PEV")) %>%
  mutate(test = ifelse(grepl("PEV", process), "uniformity", "ROI"))

data_assoc_gathered = data_assoc %>%
  gather(process, value_std, contains("ROI"), contains("PEV")) %>%
  mutate(test = ifelse(grepl("PEV", process), "association", "ROI"))

data_combined = bind_rows(data_unif_gathered, data_assoc_gathered) %>%
  unique()

head(data_combined)

large coefficient models

specify modeling approach

  • Also include the following subsets: stimulus type (food, meal, snack, dessert) and contrast (nature or rest)
setup_specs(y = c("bmi", "fat", "enact_prop"),
            x = c("reward_ROI", "reward_PEV",
                  "value_ROI", "value_PEV",
                  "cognitive_control_ROI", "cognitive_control_PEV",
                  "craving_PEV", "craving_regulation_PEV"),
            controls = c("reward_ROI", "reward_PEV",
                         "value_ROI", "value_PEV",
                         "cognitive_control_ROI", "cognitive_control_PEV",
                         "craving_PEV", "craving_regulation_PEV"),
            model = c("lm")) %>%
  filter(!x == controls)

define common variables

data = data_combined
var = c("reward_ROI", "reward_PEV",
        "value_ROI", "value_PEV",
        "cognitive_control_ROI", "cognitive_control_PEV",
        "craving_PEV", "craving_regulation_PEV")

run SCAs

BMI

# define outcome
outcome = "bmi"

# run sca
output = run_sca(data = data, var = var, outcome = outcome, boot = FALSE)

bmi_output_plot = output$results %>%
  mutate(`a priori` = ifelse(x == "value ROI" & controls == "reward ROI" &
                               stimulus == "food" & contrast == "nature", "ROI model", 
                      ifelse(x == "reward ROI" & controls == "value ROI" &
                               stimulus == "food" & contrast == "nature", "ROI model", 
                      ifelse(x == "craving regulation PEV" & controls == "no covariates" & 
                               stimulus == "food" & contrast == "nature", "PEV model", "not included"))))
  
# plot
plot_sca(data = bmi_output_plot, combined = FALSE, choices = c("x", "test", "stimulus", "contrast", "a priori"))

(bmi_compare = plot_sca_compare(data = bmi_output_plot, pointrange = FALSE))

plot_decisiontree(bmi_output_plot, legend = TRUE)

% fat

# define outcome
outcome = "fat"

# run sca
output = run_sca(data = data, var = var, outcome = outcome, boot = FALSE, n_shuffles = n_shuffles, rerun = FALSE)

fat_output_plot = output$results %>%
  mutate(`a priori` = ifelse(x == "reward ROI" & controls == "no covariates" &
                               stimulus == "food" & contrast == "nature", "ROI model", 
                      ifelse(x == "craving PEV" & controls == "no covariates" & 
                               stimulus == "food" & contrast == "nature" & test == "association", "PEV model", "not included")),
         y = ifelse(grepl("%", y), "body fat percentage", y))

# plot
plot_sca(data = fat_output_plot, combined = FALSE, choices = c("x", "test", "stimulus", "contrast", "a priori"))

(fat_compare = plot_sca_compare(data = fat_output_plot, pointrange = FALSE))

plot_decisiontree(fat_output_plot, legend = TRUE)

enactment

# define outcome
outcome = "enact_prop"

# run sca
output = run_sca(data = data, var = var, outcome = outcome, boot = FALSE, n_shuffles = n_shuffles, rerun = FALSE)

enact_output_plot = output$results %>%
  mutate(`a priori` = ifelse(x == "value ROI" & controls == "no covariates" &
                               stimulus == "food" & contrast == "nature", "ROI model", 
                      ifelse(x == "reward PEV" & controls == "no covariates" & 
                               stimulus == "food" & contrast == "nature" & test == "association", "PEV model", "not included")))

# plot
plot_sca(data = enact_output_plot, combined = FALSE, choices = c("x", "test", "stimulus", "contrast", "a priori"))

(enact_compare = plot_sca_compare(data = enact_output_plot, pointrange = FALSE))

plot_decisiontree(enact_output_plot, legend = TRUE)

combined plots

bmi_sca = plot_sca(data = bmi_output_plot, combined = FALSE, labels = c("A", "B"), title = TRUE,
                   point_size = .2, point_alpha = 1, ci_size = .1, ci_alpha = .5, text_size = 18,
                   remove_facet = TRUE,
                   choices = c("x", "test", "stimulus", "contrast", "a priori"))
fat_sca = plot_sca(data = fat_output_plot, combined = FALSE, labels = c("C", "D"), title = TRUE,
                   point_size = .2, point_alpha = 1, ci_size = .1, ci_alpha = .5, text_size = 18,
                   remove_y = TRUE, remove_facet = TRUE,
                   choices = c("x", "test", "stimulus", "contrast", "a priori"))
enact_sca = plot_sca(data = enact_output_plot, combined = FALSE, labels = c("E", "F"), title = TRUE,
                     point_size = .2, point_alpha = 1, ci_size = .1, ci_alpha = .5, text_size = 18,
                     remove_y = TRUE,
                     choices = c("x", "test", "stimulus", "contrast", "a priori"))

cowplot::plot_grid(bmi_sca, fat_sca, enact_sca, ncol = 3, rel_widths = c(.38, .31, .31))

bmi_compare = plot_sca_compare(data = bmi_output_plot, pointrange = FALSE, title = TRUE,
                               labels = c("A", "B"), rel_heights = c(.5, .5), text_size = 18,
                               n_rows = 2, 
                               sig = c("reward ROI", "value ROI", "cognitive control PEV",
                                       "craving PEV", "craving regulation PEV"))
fat_compare = plot_sca_compare(data = fat_output_plot, pointrange = FALSE, title = TRUE,
                               labels = c("C", "D"), rel_heights = c(.5, .5), text_size = 18,
                               remove_y = TRUE, n_rows = 2,
                               sig = c("value ROI", "craving PEV"))
enact_compare = plot_sca_compare(data = enact_output_plot, pointrange = FALSE, title = TRUE,
                                 labels = c("E", "F"), rel_heights = c(.5, .5), text_size = 18,
                                 remove_y = TRUE, n_rows = 2,
                                 sig = c("reward ROI", "reward PEV", "value ROI",
                                         "value PEV", "craving PEV", "cognitive control ROI"))

cowplot::plot_grid(bmi_compare, fat_compare, enact_compare, ncol = 3)

neural covariate models

define common variables

data = data_combined
n_shuffles = 500
outcome = c("bmi", "fat", "enact_prop")

run SCAs

reward ROI

# define variable
var = "reward_ROI"

# run SCA
reward_ROI_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)

reward_ROI_output$summary 
# stats
(reward_ROI = boot_stats(data = reward_ROI_output))
# plot
plot_sca(data = reward_ROI_output$results, combined = TRUE, title = TRUE)

reward PEV

# define variable
var = "reward_PEV"

# run SCA
reward_PEV_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
reward_PEV_output$summary 
# stats
(reward_PEV = boot_stats(data = reward_PEV_output))
# plot
plot_sca(data = reward_PEV_output$results, combined = TRUE, title = TRUE)

value ROI

# define variable
var = "value_ROI"

# run SCA
value_ROI_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
value_ROI_output$summary 
# stats
(value_ROI = boot_stats(data = value_ROI_output))
# plot
plot_sca(data = value_ROI_output$results, combined = TRUE, title = TRUE)

value PEV

# define variable
var = "value_PEV"

# run SCA
value_PEV_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
value_PEV_output$summary 
# stats
(value_PEV = boot_stats(data = value_PEV_output))
# plot
plot_sca(data = value_PEV_output$results, combined = TRUE, title = TRUE)

cognitive control ROI

# define variable
var = "cognitive_control_ROI"

# run SCA
cognitive_control_ROI_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
cognitive_control_ROI_output$summary 
# stats
(cognitive_control_ROI = boot_stats(data = cognitive_control_ROI_output))
# plot
plot_sca(data = cognitive_control_ROI_output$results, combined = TRUE, title = TRUE)

cognitive control PEV

# define variable
var = "cognitive_control_PEV"

# run SCA
cognitive_control_PEV_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
cognitive_control_PEV_output$summary 
# stats
(cognitive_control_PEV = boot_stats(data = cognitive_control_PEV_output))
# plot
plot_sca(data = cognitive_control_PEV_output$results, combined = TRUE, title = TRUE)

craving PEV

# define variable
var = "craving_PEV"

# run SCA
craving_PEV_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
craving_PEV_output$summary 
# stats
(craving_PEV = boot_stats(data = craving_PEV_output))
# plot
plot_sca(data = craving_PEV_output$results, combined = TRUE, title = TRUE)

craving regulation PEV

# define variable
var = "craving_regulation_PEV"

# run SCA
craving_regulation_PEV_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
craving_regulation_PEV_output$summary 
# stats
(craving_regulation_PEV = boot_stats(data = craving_regulation_PEV_output))
# plot
plot_sca(data = craving_regulation_PEV_output$results, combined = TRUE, title = TRUE)

plot combined SCAs

reward

limits = c(min(reward_PEV_output$results$conf.low), max(reward_PEV_output$results$conf.high))

a = plot_sca(data = mutate(reward_ROI_output$results, 
                           controls = ifelse(controls == "reward PEV", "reward PEV / ROI", controls)), 
             combined = TRUE, title = TRUE, labels = c("A", "B"),
             choices = c("y", "test", "stimulus", "contrast", "controls"),
             remove_facet = TRUE, text_size = 18, limits = limits)

b = plot_sca(data = reward_PEV_output$results, combined = TRUE, title = TRUE, labels = c("C", "D"),
         choices = c("y", "test", "stimulus", "contrast", "controls"),
         remove_y = TRUE, text_size = 18, limits = limits)

cowplot::plot_grid(a, b, ncol = 2, rel_widths = c(.54, .46))

value

limits = c(min(value_PEV_output$results$conf.low), max(value_PEV_output$results$conf.high))

a = plot_sca(data = mutate(value_ROI_output$results, 
                           controls = ifelse(controls == "value PEV", "value PEV / ROI", controls)), 
             combined = TRUE, title = TRUE, labels = c("A", "B"),
             choices = c("y", "test", "stimulus", "contrast", "controls"),
             remove_facet = TRUE, text_size = 18, limits = limits)

b = plot_sca(data = value_PEV_output$results, combined = TRUE, title = TRUE, labels = c("C", "D"),
         choices = c("y", "test", "stimulus", "contrast", "controls"),
         remove_y = TRUE, text_size = 18, limits = limits)

cowplot::plot_grid(a, b, ncol = 2, rel_widths = c(.53, .47))

cognitive control

limits = c(min(cognitive_control_ROI_output$results$conf.low), max(cognitive_control_PEV_output$results$conf.high))

a = plot_sca(data = mutate(cognitive_control_ROI_output$results, 
                           controls = ifelse(controls == "cognitive control PEV", "cognitive control PEV / ROI", controls)), 
             combined = TRUE, title = TRUE, labels = c("A", "B"),
             choices = c("y", "test", "stimulus", "contrast", "controls"),
             remove_facet = TRUE, text_size = 18, limits = limits)

b = plot_sca(data = cognitive_control_PEV_output$results, combined = TRUE, title = TRUE, labels = c("C", "D"),
         choices = c("y", "test", "stimulus", "contrast", "controls"),
         remove_y = TRUE, text_size = 18, limits = limits)

cowplot::plot_grid(a, b, ncol = 2, rel_widths = c(.57, .43))

craving & craving regulation

limits = c(min(craving_PEV_output$results$conf.low), max(craving_PEV_output$results$conf.high))

a = plot_sca(data = mutate(craving_PEV_output$results, 
                           controls = ifelse(controls == "craving regulation PEV", "craving regulation / craving PEV", controls)), 
             combined = TRUE, title = TRUE, labels = c("A", "B"),
             choices = c("y", "test", "stimulus", "contrast", "controls"),
             remove_facet = TRUE, text_size = 18, limits = limits)

b = plot_sca(data = craving_regulation_PEV_output$results, combined = TRUE, title = TRUE, labels = c("C", "D"),
         choices = c("y", "test", "stimulus", "contrast", "controls"),
         remove_y = TRUE, text_size = 18, limits = limits)

cowplot::plot_grid(a, b, ncol = 2, rel_widths = c(.57, .43))

combined table

bind_rows(reward_ROI, reward_PEV, value_ROI, value_PEV, cognitive_control_ROI, cognitive_control_PEV, craving_PEV, craving_regulation_PEV) %>%
  kable(format = "pandoc", digits = 2)
x BMI Mdn [95% CI] BMI Mdn p BMI Negative N BMI Negative p BMI Positive N BMI Positive p body fat percentage Mdn [95% CI] body fat percentage Mdn p body fat percentage Negative N body fat percentage Negative p body fat percentage Positive N body fat percentage Positive p enactment Mdn [95% CI] enactment Mdn p enactment Negative N enactment Negative p enactment Positive N enactment Positive p
reward ROI -0.12 [-0.14, -0.11] < .001 59 .012 0 1.000 -0.03 [-0.04, -0.02] .136 0 1.000 0 1.000 0.17 [0.15, 0.18] < .001 0 1.000 48 .014
reward PEV -0.00 [-0.01, 0.01] .392 12 .216 4 .386 -0.02 [-0.03, -0.01] .152 10 .226 0 1.000 0.20 [0.18, 0.21] < .001 4 .412 128 < .001
value ROI 0.07 [0.04, 0.08] .002 0 1.000 73 .004 0.04 [0.03, 0.05] .032 0 1.000 36 .020 0.06 [0.05, 0.08] .036 0 1.000 33 .022
value PEV 0.03 [0.02, 0.04] .038 12 .208 17 .190 0.01 [0.01, 0.02] .176 5 .424 5 .390 0.17 [0.15, 0.19] < .001 2 .432 106 < .001
cognitive control ROI 0.02 [-0.00, 0.04] .210 0 1.000 18 .106 0.07 [0.05, 0.09] < .001 0 1.000 15 .206 -0.15 [-0.19, -0.13] < .001 58 .010 0 1.000
cognitive control PEV -0.07 [-0.09, -0.06] < .001 45 .014 0 1.000 -0.02 [-0.04, -0.01] .104 8 .370 2 .474 -0.00 [-0.02, 0.00] .478 15 .212 16 .190
craving PEV 0.05 [0.04, 0.06] .006 2 .448 34 .026 0.04 [0.03, 0.05] .010 2 .420 33 .034 0.05 [0.03, 0.07] .038 14 .202 41 .018
craving regulation PEV 0.09 [0.08, 0.09] < .001 0 1.000 59 .004 0.05 [0.05, 0.06] < .001 0 1.000 8 .422 -0.00 [-0.01, 0.01] .476 6 .428 9 .326